MSIS Version Evaluation- Multiple Arcs¶
Details:¶
This notebook shows GEODYN II output for from three versions of the MSIS models for thermospheric density: MSISe86, MSISe00, MSISe2.0. The beginning of the notebook contains the analysis that can be gathered from assessing the plots that are demonstrated below.
Run Parameters:
Satellite: Starlette Observation Datatype: SLR Accel9 Adjustments: Off Density Models: NRLMSISe86, NRLMSISe00, MSIS 2.0
Data arcs in this analysis:
Arc 1 -------- 14 days ------ 03/09/14 - 03/09/28 Arc 2 -------- 14 days ------ 03/09/28 - 03/10/12 Arc 3 -------- 14 days ------ 03/10/12 - 03/10/26 Arc 4 -------- 14 days ------ 03/10/26 - 03/11/09 Arc 5 -------- 14 days ------ 03/11/09 - 03/11/23 Arc 6 -------- 14 days ------ 03/11/23 - 03/12/07 Arc 7 -------- 14 days ------ 03/12/07 - 03/12/21 Arc 8 -------- 14 days ------ 03/12/21 - 04/01/5
Analysis¶
Background:¶
Starlette Orbit Info:
APOGEE = 1110727.924 M
PERIGEE= 810723.637 M
ECCEN. = 0.0204
INCLIN = 49.8080 DEG
The slight eccentricity accounts for much of the short time scale variation that can be seen in the density (Seen here: Density Plot).
Doing an orbit average (Seen here: Orbit Avg’d Plot) removes much of the variation due to eccentricity such that the remaining variation is caused primarily by geomagnetic and solar flux variation.
Increased density from Geomagnetic and solar activity :
the Orbit Avg’d Plot) also shows the geomagnetic and solar activity in the form of the Ap and F10.7, respectively.
Because the satellite is in a mid-latitude inclination orbit, the spike in density seen in Arcs 4 and 5 is could be more due the increase in solar flux than due to increased Ap.
Relative model performance:¶
We have to look at a combination of different factors. This notebook includes the following: 1. RMS of fit on final iteration (want small) 2. Cd corrections from a priori (want small) - In general, the drag coefficient will be adjusted (corrected by some amount) to account for variations seen in the mass density provided by the different models. This makes the Cd (in this run scenario) a physically meaningful parameter to assess the POD. 3. Observation residuals (want to be small, but don’t rely too heavily on this)
*The mean in the residuals (or the mean by station) is not important, it is very small in these runs*
Adjusted Cd (apriori value = 2.2) - At this altitude regime (and over this epoch), the MSISe86 model has the smallest corrections from a priori at a mean percent difference of 18.4%. MSISe00 has slightly higher corrections at 19.1 %. MSISe2.0 has a mean percent difference from a priori of 20.25%. - In stormtime scenarios, the Cd respond by correcting the mean density provided by the different models. - At this altitude regime, the MSIS2.0 has the largest correction during storm time events. The end of Arc 5 shows a particularly large correction in which all models correct the Cd by more than 200% (m2=~280%, m00=~240%, m86=~205%). - from a POD perspective, the Cd does a good job of accounting for the density variations. Its not clear to me at this stage what the physical significance of the corrections is… [ASK QUESTION]
RMS of Fit
MSISe86 also appears to have the smallest RMS of fit of these model versions. RMS = .0922857
MSISe00 has a slighly higher RMS at .0923714
MSISe86 has the highest RMS of these models at .0931429
It is worth noting that since this is at such a high altitude, the many improvements that were made to the MSIS models in version 00 and 2.0 may not be impacting these results in a POD sense. The data that went into MSIS2.0 was mostly from satellites and data sets derived at lower altitudes. We will need to repeat the assessment methods demonstrated here for a lower altitude satellite geodyn run and compare the results.
Load the simulations over all arcs¶
[1]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px
import copy
import sys
import numpy as np
[9]:
%load_ext autoreload
%autoreload 2
import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/')
from PYGEODYN import Pygeodyn
#------ A dictionary containing the run parameters ------
run_params1 = {}
run_params1['arc'] = [ '030914_2wk',
'030928_2wk',
'031012_2wk',
'031026_2wk',
'031109_2wk',
'031123_2wk',
'031207_2wk'
]
run_params1['satellite'] = 'starlette'
run_params1['den_model'] = 'msis2'
run_params1['SpecialRun_name'] = ''
run_params1['verbose'] = False
run_params1['action'] = 'read'
Obj_Geodyn1 = Pygeodyn(run_params1)
Obj_Geodyn1.getData()
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
=================================================================
Initializing PYGEODYN to ...
... READ GEODYN output
=================================================================
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
<ipython-input-9-12b036d5271d> in <module>
24
25 Obj_Geodyn1 = Pygeodyn(run_params1)
---> 26 Obj_Geodyn1.getData()
27
/data/geodyn_proj/pygeodyn/utils_pygeodyn_develop/PYGEODYN_Read.py in getData(self)
2003 os.chdir(self.path_to_model+'DENSITY/')
2004 os.system('bunzip2 -v '+'*')
-> 2005 os.chdir(self.path_to_model+'ORBITS/')
2006 os.system('bunzip2 -v '+'*')
2007
FileNotFoundError: [Errno 2] No such file or directory: '/data/data_geodyn/results/icesat2/msis2/msis2_acceloff/ORBITS/'
[3]:
params_msis2 = copy.deepcopy(run_params)
params_msis2['den_model'] = 'msis2'
obj_msis2 = pyGeodyn_Readers( params_msis2 )
obj_msis2.getData_All()
Calling pygeodyn with multiple arcs...
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031207_2wk.goco05s
[4]:
params_msis00 = copy.deepcopy(run_params)
params_msis00['den_model'] = 'msis00'
obj_msis00 = pyGeodyn_Readers( params_msis00 )
obj_msis00.getData_All()
Calling pygeodyn with multiple arcs...
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031207_2wk.goco05s
[5]:
params_msis86 = copy.deepcopy(run_params)
params_msis86['den_model'] = 'msis86'
obj_msis86 = pyGeodyn_Readers( params_msis86 )
obj_msis86.getData_All()
Calling pygeodyn with multiple arcs...
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis86/msis86_acceloff/.../st031207_2wk.goco05s
Set the plot configurations¶
[6]:
config = dict({
'displayModeBar': True,
'responsive': False,
# 'staticPlot': True,
'displaylogo': False,
'showTips': False,
})
Coefficient of Drag Adjustments¶
Details and Background:¶
the C_d is an adjusted parameter in GEODYN. Each iteration, it is estimated and adjusted to improve the OD simulation
C_d is treated as a constant in GEODYN (for any iteration). The factor C_d varies slightly with satellite shape and atmospheric composition. However, for any geodetically useful satellite, it may be treated as a satellite dependent constant. If the Cd is time dependent, then it is constant over the designated time intervals.
Cd Analysis¶
A few key features can be seen in the above plot: 1. When there is some sort of storm event early in the arc, the Cd responds by correcting the mean density provided by the different models. 2. The MSISe00 corrections are less at this set of altitudes as opposed to those of m86 and m2.0. 3. Overall, though from a POD perspective, they seem to do a good job accounting for density variations. - Frank states that the onus is to show there is some real storm event going on over these several days.
Adjusted Cd Plot (timeseries)¶
[7]:
%load_ext autoreload
%autoreload 2
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/utils_pygeodyn')
from pygeodyn_Analysis import plot_cd_and_percdiff_from_apriori
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
[8]:
fig = make_subplots(
rows=2, cols=1,
subplot_titles=(["Timeseries of Cd", "Percent difference from a priori (Cd=2.2)"]),
vertical_spacing = 0.1,
)
fig = plot_cd_and_percdiff_from_apriori(fig, obj_msis2, 0)
fig = plot_cd_and_percdiff_from_apriori(fig, obj_msis00, 1)
fig = plot_cd_and_percdiff_from_apriori(fig, obj_msis86, 2)
iplot(fig)
fig.show(config=config)
Scale Density with the Cd scale factor:¶
The percent difference of the estimated/adjusted Cd from the a priori Cd can be seen as a scale factor that is being applied to the drag equation to account for unknown or poorly modeled variations in the neutral density. By multiplying the density by the scale factor outputed by the percent change in the Cd from apriori, we get a good estimation of the “more true” density.
[9]:
%load_ext autoreload
%autoreload 2
from pygeodyn_Analysis import plot_ScaleDensity_with_CdScaleFactor
fig = make_subplots(rows=2, cols=1,
subplot_titles=(["Regular and Scaled Density", "Scale Factor"]),
vertical_spacing = 0.1,
specs=[
[{"secondary_y": False}],
[{"secondary_y": True}], ])
fig = plot_ScaleDensity_with_CdScaleFactor(fig, obj_msis2, 0, 200)
# iplot(fig)
fig.show(config=config)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Density Output¶
The density output is sourced from the
DRAG.f90subroutine.
Density output Analysis¶
Starlette is in somewhat of an eccentric orbit, roughly 800 km x 1100 km orbit - so that accounts for the change in density along the orbit.
while this is helpful to see, one really needs to demonstrate these results over multiple arcs.
The relative densities are different at Apogee and perigee.
Apogee : MSISe86 and MSISe00 have the higher densities
Perigee : MSIS2.0 has the higher density
Density timeseries plot¶
[10]:
%load_ext autoreload
%autoreload 2
from pygeodyn_Analysis import plot_density_decimated
fig = make_subplots(
rows=1, cols=1,
subplot_titles=(['Every 200th point']),
vertical_spacing = 0.1,
)
fig = plot_density_decimated(fig, obj_msis2, 0, 200)
fig = plot_density_decimated(fig, obj_msis00, 1, 200)
fig = plot_density_decimated(fig, obj_msis86, 2, 200)
# iplot(fig)
fig.show(config=config)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Orbit averaged Density Plot¶
[11]:
# %load_ext autoreload
# %autoreload 2
# from pygeodyn_Analysis import plot_density_orbit_avg_with_fluxes
# fig = make_subplots(rows=3, cols=1,
# subplot_titles=(['Orbit Averaged Density', 'Ap', 'F10.7']),
# horizontal_spacing = 0.05,
# vertical_spacing = 0.1,
# )
# fig = plot_density_orbit_avg_with_fluxes(fig, obj_msis2, 0)
# fig = plot_density_orbit_avg_with_fluxes(fig, obj_msis00, 1)
# fig = plot_density_orbit_avg_with_fluxes(fig, obj_msis86, 2)
## iplot(fig)
# fig.show(config=config)
Composite Plots, density, cd, activity indicies¶
[12]:
%load_ext autoreload
%autoreload 2
from pygeodyn_Analysis import plot_composite_density_cd_fluxes
fig = make_subplots(rows=3, cols=1,
subplot_titles=(['Orbit Averaged Density','Adjusted Cd', 'Geomagnetic and Solar Activity']),
horizontal_spacing = 0.05,
vertical_spacing = 0.1,
specs=[[{"secondary_y": False}],
[{"secondary_y": False}],
[{"secondary_y": True}], ])
fig = plot_composite_density_cd_fluxes(fig, obj_msis2, 0)
fig = plot_composite_density_cd_fluxes(fig, obj_msis00, 1)
fig = plot_composite_density_cd_fluxes(fig, obj_msis86, 2)
# iplot(fig)
fig.show(config=config)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Residuals – obervations¶
Details:¶
the residuals being output by geodyn should be rescaled to CM instead of METER. Most of your stations have cm level (or better) data precision.
The mean in the residuals (or the mean by station) is not important - it is very small in these runs - unless you turn out to have a station with a large systematic range bias - due to some temporary anomaly.
If you look at the setups, you’ll see sets of “MBIAS” cards – these adjust different sets of range biases by station per arc or by station per data pass based on recommendations from the ILRS. They are adjusted as a normal part of the POD process.
Analysis¶
You see the bowtie effect in the residuals above. That is a standard effect of least squares - where fits seem the best in the center of the arc.
Observed Residuals timeseries plot¶
[13]:
%load_ext autoreload
%autoreload 2
from pygeodyn_Analysis import plot_residuals_observed
fig = make_subplots(rows=1, cols=1,
# subplot_titles=(['Orbit Averaged Density', 'Ap', 'F10.7']),
# horizontal_spacing = 0.05,
# vertical_spacing = 0.1,
)
fig = plot_residuals_observed(fig, obj_msis2, 0)
fig = plot_residuals_observed(fig, obj_msis00, 1)
fig = plot_residuals_observed(fig, obj_msis86, 2)
# iplot(fig)
fig.show(config=config)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Residual Measurement Summary plots (RMS of fit)¶
[14]:
%load_ext autoreload
%autoreload 2
from pygeodyn_Analysis import plot_residual_meas_summary
fig = make_subplots(rows=2, cols=1,
subplot_titles=(["Mean Residuals", 'RMS of Fit']),
vertical_spacing = 0.1,
)
fig = plot_residual_meas_summary(fig, obj_msis2, 0)
fig = plot_residual_meas_summary(fig, obj_msis00, 1)
fig = plot_residual_meas_summary(fig, obj_msis86, 2)
# iplot(fig)
fig.show(config=config)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
msis2 RMS 9.31429e-02
msis2 Resid -3.71429e-02
msis00 RMS 9.23714e-02
msis00 Resid -3.28571e-02
msis86 RMS 9.22857e-02
msis86 Resid -3.28571e-02
Residuals – station summary¶
I have been told that these are not helpful for assessing density models.
unless to identify a really bad station (they can crop up from time to time because of a station anomaly). If one finds a bad station, then one can just delete the data. In this example, it looks like Riyadh (RIYA7832) seems to have a mean bias of 2 cm?. That is a bit large and might be a sign of station coordinate error or some other station-dependent error. One could solve for an MBIAS for that station over the arc to take care of that small problem - you would have to check that an MBIAS for that station is not already adjusted.
TODO: NEED TO ASK FL WHAT EXACTLY THESE DO